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Abstract. Self-sustained dynamical phases of living matter can exhibit remarkable 
similarities over a wide range of scales, from mesoscopic vortex structures in microbial 
suspensions and motility assays of biopolymers to turbulent large-scale instabilities in 
flocks of birds or schools offish. Here, we argue that, in many cases, the phenomenology 
of such active states can be efficiently described in terms of fourth- and higher-order 
partial differential equations. Structural transitions in these models can be interpreted 
as Landau-type kinematic transitions in Fourier (wavenumber) space, suggesting 
that microscopically different biological systems can share universal long-wavelength 
features. This general idea is illustrated through numerical simulations for two classes 
of continuum models for incompressible active fluids: a Swift-Hohenberg-type scalar 
field theory, and a minimal vector model that extends the classical Toner- Tu theory and 
appears to be a promising candidate for the quantitive description of dense bacterial 
suspensions. We also discuss briefly how microscopic symmetry-breaking mechanisms 
can enter macroscopic continuum descriptions of collective microbial motion near 
surfaces. 
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1. Introduction 

Simple and complex life forms can exhibit remarkably similar collective behaviors over 
a wide range of length and time scales [U El E]- Well-known examples are flocking 
phenomena in swarms of birds jlj and self-sustained turbulent phases in schools of 
fish |5] that share several qualitative features with the meso-scale dynamics in bacterial 
suspensions [6j [7J EJ [9] and films [10J [III [12] • Over the past two decades, intense 
efforts have been made to understand the phenomenology of microbial and other 
active fluids, but in spite of substantial progress it is still not entirely clear which 
of their characteristics are universal or system-specific [21 [321 [33], and which classes 
of dynamical equations are capable of providing adequate minimal descriptions. In 
recent years, a considerable number of continuum models for active systems have been 
proposed [21 [151 [TBI fT7| fTH^ l2T]. but most of them have yet to be tested against 

experiments [221 [23]. Many of those theories focus on the couplings between two or more 
order-parameters (concentration, solvent velocity, orientation fields, etc.) and typically 
involve a large number of parameters, thus making comparison with experimental data 
very difficult. To exploit current and future progress in experimental imaging and 
tracking techniques j5j EU [251 EEJ E7J , and to understand better the general ordering 
principles that govern active matter [H [2], it will be necessary to identify tractable 
minimal models that not only capture the essential instability mechanisms but also 
allow for quantitative comparison with experiments. 

In this paper, we will analyze two such minimal continuum theories for active 
suspensions by focussing on generic structural properties and stressing formal analogies 
with classical phase transitions. Our approach is based on the hypothesis that dynamical 
transitions in many internally or externally driven systems, such as microbial [HI [TJ 
El ESI ESI EH] or vibrated colloidal suspensions [291 E0], can be phenomenologically 
modeled as Landau-type transitions in Fourier (wave-number) space, which suggests 
that minimal hydrodynamic descriptions of active matter can be obtained in terms of 
higher-than-second-order partial differential equations (PDEs). Higher-order PDEs have 
been previously derived and studied for a wide range of nonlinear structure formation 
phenomena [221 EH E21 ESI El], including Rayleigh-Benard convection [33], polymer 
and vesicle dynamics [36], quasi-crystal formation [37] and theories of ionic liquids [38J. 
However, to our knowledge, models of this type have rarely been considered in the 
context of microbial suspensions [39J. Therefore, one of our main objectives here is 
to draw attention to the possibility that one can obtain useful, testable continuum 
theories of bacterial and other active fluids, by restricting field variables to a minimal 
set of experimentally accessible order-parameters but admitting fourth- and higher-order 
spatial derivatives. 

The basic idea is readily summarized as follows: It is well-known that the 
incompressible Navier-Stokes equation is capable of describing the dissipative flow 
dynamics v(t, x) of a wide range of conventional 'passive' fluids, regardless of their exact 
microscopic composition. A main reason for this is that these systems behave similarly 
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at long- wavelengths (small- wave numbers), so that the leading order viscous dissipation 
can be described by a term T Av in the field equations, with partial information 
about the microstructure being retained in the viscosity coefficient To (throughout, 
A = V 2 denotes the Laplacian). Transforming to Fourier-space, the dissipative term 
yields a simple quadratic 'dispersion' relation ~ r |k| 2 , which represents the dominant 
contribution in a systematic small- wavenumber expansion and leads to damping in the 
absence of external stimuli. By contrast, in active fluids, viscous dissipation competes 
with internal or external energy input and, in principle, one cannot exclude that higher- 
order contributions of the form r |k| 2 + r 2 |k| 4 + . . . become relevant as well. In fact, 
they will certainly be needed to ensure stability if, due to the complex interplay of 
nonlinear interactions and energy input, the coefficient To should change its sign. The 
inclusion of higher-order terms in the Fourier-space expansions is analogous to the well- 
known Landau-expansion of order-parameter potentials and, accordingly, sign-changes 
in the coefficients T n can give rise to Landau-type kinematic phase transitions. When 
going back to position space, terms |k| 4 , |k| 6 , . . . will transform into higher-order spatial 
derivatives A 2 , A 3 , .... The inclusion of such terms makes the theory successively more 
non-local. The restriction to functions of |k| is dictated by isotropy; in principle, one 
could also study odd and fractional powers of |k| but this would go beyond the scope of 
the present paper. 

From the preceding considerations, it seems plausible that a systematic 
characterization of active fluids in terms of their asymptotic small wave-number 
expansions can help to distinguish specific from universal properties, thereby providing 
a basis for more systematic classification schemes similar to those for thermodynamic 
equilibrium phases in classical fluids or spin systems. Moreover, this analogy- 
driven approach promises analytically tractable models of active suspensions that are 
considerably simpler than many of the currently studied (potentially more accurate) 
multi-component theories and will hopefully enable quantitative comparisons with 
experiments in the near future. In the present paper, we shall focus on theoretical aspects 
of fourth-order continuum models, starting with the simplest case, which is given by a 
Swift-Hohenberg-type scalar or pseudo-scalar field theory [311 EH] • This model is used 
as a basic example to illustrate how microscopic symmetry-breaking mechanisms [ID] 
can enter macroscopic continuum descriptions of microbial motion near surfaces [28J. 
Subsequently, we will generalize to non-scalar order-parameters by considering a minimal 
vector theory for incompressible active suspensions. The resulting flow model extends 
the seminal Toner- Tu theory [131 HB] and appears to be a promising candidate for the 
quantitive description of highly concentrated bacterial fluids [23] ■ In the subsequent 
discussion, we use results from two-dimensional (2D) continuum simulations to illustrate 
selected dynamical properties of the different models in more detail. 
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2. (Pseudo) scalar order-parameter theory 

The minimal model considered in this section belongs to the class of generalized 
Swift-Hohenberg theories [29| 135] . Our motivation for prepending a brief discussion 
of this well-known model here is two-fold: It is helpful to recall some of its basic 
properties before considering the generalization to vectorial order-parameters. This 
model is also useful for illustrating how microscopic symmetry-breaking mechanisms [40J 
can be incorporated into macroscopic descriptions of experimentally relevant microbial 



systems [28J, as discussed in Section 2.4 below. 



2.1. Model equations 

We consider the simplest isotropic fourth-order model for a non-conserved scalar or 
pseudo-scalar order-parameter if)(t,x), given by 

a^ = F(V)+7oAV- 72 A 2 V, (1) 

where dt = d/dt denotes the time derivative, and A = V 2 is the (^-dimensional 
Laplacian. The force F is derived from a Landau-potental U(i/j) 

F = -^, U W = a -^ + b -^+ C ^\ (2) 

and the derivative terms on the rhs. of Equation ([!]) can also be obtained by variational 
methods from a suitably defined energy functional. In the context of active suspensions, 
■0 could, for example, quantify local energy fluctuations, local alignment, phase 
differences, or vorticity. We will assume throughout that the system is confined to 
a finite spatial domain Q C M d of volume 

= / d d x, (3) 
Jn 

adopting with periodic boundary conditions in simulations. 

For completeness, one should note that in the case of a conserved order-parameter 
field q the field equations would either have to take the current-form d t g = — V -3(g) 
or, alternatively, one could implement conservation laws globally by means of Lagrange 
multipliers [31j] . For example, for a dynamics similar to that of Equation ([I]) and a 
simple global 'mass' constraint 



M = / d x q = const, 
Jn 

the Lagrange-multiplier approach yields the non-local equations of motions 
d t Q = F(g) + 7oA£ - 7 2 A 2 £ - Ai, 
Al = W\ I ^ [ F (£)+7oA£- 72 A 2 £] . 

I | */ $6 

In the remainder of this section, however, we shall focus on the local dynamics defined 
by Equations ([T]) and ([2]), since this well-known example will be a useful reference point 
for the discussion of the vector model in Section [3] 
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2.2. Linear stability 

The fixed points of Equation Q are determined by the zeros of the force F(ip), 
corresponding to the minima of the potential U, yielding ipo = and 



b b 2 a „ , 9 

^ = -_±^_-_, if ^>4ac. (4) 

Linearization of Equation near ^ f° r small perturbations ip — eo ex P( — °"o^ — «k • x) 
gives 

<7 (k) = a + 7o |k| 2 + 72 |k| 4 . (5) 

Similarly, one finds for if) = ip± + e± exp(— o~±t — ik ■ x) 

<j±(k) = -(2a + 6V±) + 7o|k| 2 + 7 2 |k| 4 . (6) 

The unusual sign-convention in the exponential of the perturbation ansatz was so chosen 
as to emphasize the formal similarity of Equations ^ and ^ with the quartic Landau 
potential (j2|, i.e., modes with a < are unstable. 

From Equations ^ and (|6]), we see immediately that 72 > is required to ensure 
small-wavelength stability of the theory and, furthermore, that non-trivial dynamics 
can be expected if a and/or 70 take negative values. In particular, all three fixed points 
can become simultaneously unstable if 70 < 0. The analogy with classical Landau- 
transitions is evident if we compare ^ and ^ with the order-parameter potential U in 
Equation ^ for the symmetric case 6 = 0: Changing the sign of 70 induces a dynamical 
transition (in Fourier space), which is formally similar to the standard 'configurational' 
second-order transition j2] in the vicinity of a = 0. 



2.3. Numerical results in 2D 

We briefly illustrate the 7 -induced changes in the dynamics of the (pseudo-) scalar 
field ip(t,x.) through 2D numerical results. The discussion in this part merely serves as 



a reminder before considering symmetry-breaking in Section 2.4 



Algorithm To simulate Equations ([T]) and (|2]) in two space dimensions, we implemented 
a pseudospectral algorithm with periodic boundary conditions as commonly used in 
computational fluid dynamics [41J. The model equations are projected onto a Fourier 
space basis, and the remaining ordinary differential equations are solved numerically 
by an operator splitting method that computes the linear operator exactly [32] • The 
nonlinear terms were evaluated by applying the '2/3-rule' to suppress aliasing errors [43J. 
We simulated the model dynamics on 2D cubic grids with sizes ranging from 64 x 64 to 
256 x 256 lattice points. The solver was written in Matlab, and its numerical stability 
was verified for a wide range of parameters and space-time discretizations. Rescaled 
dimensionless variables and parameters as adopted in the simulations are summarized 
in Table [T] The rescaled time steps were typically of the order of At = 1CT 1 . All 
simulations were initiated with isotropic, randomly chosen order-parameter values. 
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model parameter 


rescaled dimensionless parameter 


a 


a t u 


b 


b 4> u U 


c 


1 


7o 


7o t u L- 2 


72 


1 



Table 1. The right column shows the rescaled dimensionless parameters used in the 
simulations of the (pseudo)scalar model from Equations ([I]) and The unit time 
is defined by the damping time-scale t u = L 4 /^, where L is the length of the 2D 
simulation box, and the order- parameter is measured in units of ^ u = l/y/t u c. 



Structural transitions Results from the numerical simulations for the order-parameter 
field ip(t,x) and two qualitatively different potentials U(ip) are summarized in Figure [I] 
In these simulation, the parameter 70 was varied between successive runs while keeping 
all other parameters fixed. To quantify changes in the quasi-stationary dynamics 
of ip as a function of 70, we measured the space-time averaged standard deviation 
u % = i^ 2 ) ~ W 2 (Figure [l^,d). Regions with = correspond to disordered 
structureless stationary states, whereas a\ > indicates the emergence of stationary 
or quasi- stationary dynamical structures. Singular points in the curve a 2 (70) signal 
qualitative changes in the order-parameter dynamics. 

In the case of a mono-stable potential (a > 0), the quantifier ex 2 undergoes a series 
of continuous transitions as 70 is lowered to negative values, see Figure [l^,. Each of 
those transitions corresponds to an increase in the number of 'stripes' that are found 
to persist for long periods of time in the simulations (Figure |lp,c). By contrast, in the 
case of bi-stable potentials (a < 0), the onset of pattern formation carries the signature 
of a first-order transition reflected by a sudden jump in cr| (Figure [ljl). However, 
while such singularities in share some formal similarities with macroscopic phase 
transitions, one could also argue that they merely signal a change in the typical number 
of excitable modes in the system. In fact, by viewing such elementary excitations as 
'quasi-particles', the structural transitions in Figure [l^,d appear to be more closely 
related to finite-systems singular points [1U HSJ H6] . 

An estimate of the critical absolute value 70 for the first 'disorder-structure' tran- 
sition can be obtained by dimensional analysis, or by equating the last two terms in 
Equations ([5j or (J6|, yielding 7q ps — (27r) 2 7 2 /L 2 . For 70 7q, increasingly more com- 
plex quasi-stationary patterns may arise (Figure [l]3,d). Structures similar to those in 
Figure [l] have been observed and widely studied [33] in granular media [29] and chemical 
systems [17] ■ In the next section, we shall demonstrate that certain aspects of collective 
microbial motion can be described within the same class of fourth-order PDEs. 
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(a) 



— i — i — i — i — i — 

mono-stable potential 
a = 0.2, 6 = 



(d) 2 




bi-stable potential 
a = -0.2, 6 = 




0.5 




Figure 1. Numerical illustration of structural transitions in the order-parameter 
ip for (a-c) mono-stable and (d-f) bi-stable potentials. (a,d) Symbols show the 
results of simulations for the first two 70-induced transitions, and lines are linear 
interpolations. Quasi-stationary space-time averages ( • ) were computed over 3000 
successive simulation time-steps (At = 0.1) after an initial relaxation period of 200 
characteristic time units t u = L 4 /^- (b,c) Snapshots of the order-parameter field ip 
at t = 500, scaled by the maximum value ip m , for a mono-stable potential U(ip) and 
homogeneous random initial conditions. After the first transition two stripes appear, 
and the number of stripes increases with the number of transitions. (e,f) Snapshots 
of the order-parameter at t = 500 for a bi-stable potential. For 70 -C — (27r) 2 72/I/ 2 , 
increasingly more complex quasi-stationary structures arise; see References |29l 147] for 
similar patterns in excited granular media and chemical reaction systems. 
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2.4- Symmetry breaking 

With regard to microbial suspensions, the minimal model ([I]) is useful for illustrating 
how microscopic symmetry-breaking mechanisms that affect the motion of individual 
organisms or cells jJOj HHJ H9J |50] can be implemented into macroscopic field equations. 
To demonstrate this, we interpret ip as a 2D pseudo-scalar vorticity fielc^jj] 

i|)eu = VAv = eijdiVj, (7) 

which is assumed to describe the flow dynamics v of a dense microbial suspension 
confined to a thin quasi-2D layer of fluid. If the confinement mechanism is top-bottom 
symmetric, as for example in a thin free-standing bacterial film [10], then one would 
expect that vortices of either handedness are equally likely. In this case, Equation (jl]) 
must be invariant under oj — > —ou, implying that U(uS) = U(—u) and, therefore, 6 = 
in Equation ([2]). Intuitively, the transformation u — > -co corresponds to a reflection of 
the observer position at the midplane of the film (watching the 2D layer from above vs. 
watching it from below). 

The situation can be rather different, however, if we consider the dynamics of 
microorganisms close to a liquid-solid interface, such as the motion of bacteria or sperms 
cells in the vicinity of a glass slide (Figure [2]). In this case, it is known that the 
trajectory of a swimming cell can exhibit a preferred handedness pH)} H8l H9j |50] . For 
example, the bacteria Escherichia coli [10] and Caulobacter [IB] have been observed 

|| 6ij denotes the Cartesian components of the Levi-Civita tensor, <9j = d/dxi for i = 1, 2, and we use 
a summation convention for equal indices throughout. 




(a) ip/tl)-_ 




0.5 1 

x/L 



Figure 2. Effect of symmetry breaking, (a) Stationary hexagonal lattice of the pseudo- 
scalar vorticity order-parameter ip = w, scaled by the maximum value tp m = L) m , 
as obtained in simulations of Equations ([T]) and ^ with b > 0, corresponding to a 
broken reflection symmetry u) — w. Blue regions correspond to clockwise motions, 
(b) Hexagonal vortex lattice formed spermatozoa of sea urchins (Strongylocentrotus 
droebachiensis) near a glass surface; from |28j adapted and reprinted with permission 
from AAAS. At high densities, the spermatozoa assemble into vortices that rotate in 
clockwise direction (inset) when viewed from the bulk fluid. 
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to swim in circles when confined near to a solid surface. More precisely, due to an 
intrinsic chirality in their swimming apparatus, these organisms move on circular orbits 
in clockwise (anticlockwise) direction when viewed from inside the bulk fluid (glass 
surface). Hence, the presence of the near-by no-slip boundary breaks the reflection 
symmetry, u -/> —u. The simplest way of accounting for this in a macroscopic continuum 
model, is to adapt the potential U(u) by permitting values b ^ in Equation The 
result of a simulation with b > is shown in Figure [2^i. In contrast to the symmtric 
case 6 = (compare Figure [lp), an asymmetric potential favors the formation of 
stable hexagonal vorticity patterns (Figure [2^,) - such self-assembled hexagonal vortex 
lattices have indeed been observed experimentally [2E] for spermatozoa of sea urchins 
(Strongylocentrotus droebachiensis) near a glass surface (Figure [2b). 



3. Vector model for an incompressible active fluid 



We now generalize the preceding considerations to identify a minimal vector-field model 
for dense microbial suspensions. Previously developed continuum theories [2J EEl HEJ [TFJ 
IT8j IT9~| l20~t |2T] of microbial fluids typically distinguish solvent concentration, bacterial 
density, solvent velocity, bacterial velocity, and various orientational order-parameter 
fields (polarization, Q-tensors, etc.). Aiming to identify a minimal hydrodynamic model, 
we construct a simplified higher-order theory by focussing exclusively on the dynamics 
of the mean bacterial velocity field v(t, x) and restricting ourselves to the incompressible 
limit. By construction, the resulting v-only theory, which is essentially a minimal Swift- 
Hohenberg-type [22] extension of the Toner- Tu model [13 [T6] , may not be applicable 
to swarming or flocking regimes, where density fluctuations are dominant, but it can 
provide a useful basis for quantitative comparisons with experiments and simulations on 
highly concentrated active suspensions [23]. In practice, v can be determined applying 
suitable coarse-graining procedures (PIV algorithms, local averaging, etc.) to discrete 
experimental or numerical velocity data [23l I5T] . 



3.1. Model equations 
Postulating incompressibility 

V • v = d t Vi = 0, (8) 
we assume that the dynamics of v is governed by the generalized Navier-Stokes equation 

(<9 t + v • V)v = -Vp -(A + C|v| 2 )v + V • E. (9) 

The pressure p(t, x) is the Lagrange multiplier for the incompressibility constraint. 
Similar to the scalar case, Equation ^ above, the (A, C)-terms in Equation ^ represent 
a quartic Landau velocity potential [21 HS1 EES] 

t/(v) = ^|v| 2 + ^|v| 4 . (10) 
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For A > and C > 0, the potential is mono-stable and the fluid is damped to a 
disordered state with v = 0. By contrast, for A < 0, Equation (10) describes a d- 
dimensional mexican-hat (sombrero) potential that favors global polar order. However, 
the fact that homogeneous polar ordering is not observed in suspensions of swimming 
bacteria J7] suggests that other instability mechanisms prevail [2D] . Mathematically, this 
means that one must either introduce additional order parameters [21 [15l [16] or make the 
theory more non-local by including higher-order derivatives via the stress-tensor [35J. 
Adopting the latter approach, we postulate that the components of the symmetric and 
traceless rate-of-strain E tensor are given by 

= T (diVj + djVi) - T 2 A (diVj + djVi) + S q^ (11) 

where 

qij = ViVj - —\\\ (12) 

is a d x rf-dimensional mean-field approximation to the Q-tensor, representing active 
stresses due to swimming (5^ is the Kronecker tensor). Although the S-term does 
not affect the linear stability of the model, general hydrodynamic considerations [T9] 
suggest that S < for pusher-swimmers like E. coli (52] or B. subtilis, whereas S > for 
puller- type microswimmers such as Chlamydomonas algae [53]. Inserting Equations (11) 
and (12) into Equation (J9|, and defining 

A = 1 - S, Ai = -S/d, (13) 

we obtain 

(d t + A v • V)v = -Vp + AiVv 2 - (A + C|v| 2 )v + r Av - T 2 A 2 v. (14) 



For T 2 = and T > 0, Equation (14) reduces to an incompressible version of 
the classical Toner- Tu model [21 [15], HE]. It is, however, the combination of the 
two T-terms with the non-variational convective derivative that turns out to be 
crucial for the formation of self-sustained quasi-chaotic flow patterns. The linear la- 
terals are reminiscent of the higher-order spatial derivatives in the classical Swift- 
Hohenberg theory [35], and Equation (14) yields a simple - if not the simplest - 
generic continuum description of turbulent meso-scale instabilities observed in dense 
bacterial suspensions [23] • More generally, Equation (14) can provide a satisfactory 
phenomenological model whenever interaction terms in more complex field theories, 
that lead to instabilities in the v-field, can be effectively approximated by a fourth- 
order Taylor expansion in Fourier space. This is likely to be the case for a wide range of 
active systems. Phrased differently, the last two terms in Equation (14) may be regarded 
as the Fourier-space analogue of the Toner- Tu driving terms, which correspond to a 
series expansion in terms of the order-parameter. Hence, similar to the higher-order 
gradient terms in the scalar theory from Equation 0, the (r , r 2 )-terms in Equation 
(14) describe intermediate-range interactions, and their role in Fourier-space is similar 
to that of the Landau potential in velocity space. 
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3.2. Linear Stability Analysis 

To support the qualitative statements in the preceding paragraph, we now perform a 
stability analysis for the 2D case relevant to the simulations discussed below, assuming 
r < and C > 0, T 2 > 0. 



The fixed points of Equations l\8ti and ( 14 ) are given by the extrema of the quartic 



velocity potential U(v). For arbitrary values of A, Equations ([8]) and (14) have a fixed 
point that corresponds to a disordered isotropic state (v,p) = (0,p ) where p is a 
constant pressure. For A < 0, an additional class of fixed points arises, corresponding 
to a manifold of globally ordered polar states (v,p) = (y ,p ), where v is constant 
vector with arbitrary orientation and fixed swimming speed |vq| = yf—AjC =: vq. 



Linearizing Equations (|8j) and (14) for small velocity and pressure perturbations 
around the isotropic state, v = e and p = p + r\ with \r\\ <C \po\, an d considering 
perturbations of the form 



0?,e) = 07,e)exp(-<7 O £-zk-x), (15) 



we find 



=ke, (16) 
a e= -^k+(A + r |k| 2 + r 2 |k| 4 )e. (17) 

Multiplying the second equation by k and using the incompressibility condition implies 
that fj = and, therefore, 

<7 (k) = A + r |k| 2 + r 2 |k| 4 . (is) 

Assuming r < and T 2 > 0, and provided that 4A < |r | 2 /r 2 , we find an unstable 
band of modes with o"o(k) < for k 2 _ < |k| 2 < k+, where 




For A < the isotropic state is generally unstable with respect to long- wavelength (i.e., 
small- 1 k|) perturbations. 

We next perform a similar analysis for the polar state (vo,Po)j which is energetically 
preferred for A < and corresponds to all active particles swimming in the same 
direction ('global order'). In this case, when considering small deviations 

v = v + e, p = Po + r], (20) 

it is useful to distinguish perturbations perpendicular and parallel to v , by writing 
e = e|| + e_i_ where v • e± = and v • €\\ = vq€\\. Without loss of generality, we may 
choose v to point along the x-axis, v = Voe x . Adopting this convention, we have 
6[| = (e[[,0) and e± = (0, e±), and to leading order 

|v| 2 ~t; 2 + 2 Wo e||. (21) 

Linearization for exponential perturbations of the form 

(V, e lh = (V, en, ej_) exp(-crt - zk • x) (22) 



Minimal continuum models of active fluids 



12 



yields 



where 



= k • e, 

a e = — i(fj — 2uoAieii)k — Me, 



m=i 2 Q A ° ] -(r |k| 2 + r 2 |k| 4 -a M )i 



(23) 
(24) 

(25) 



with I = (5ij) denoting the identity matrix. Multiplying Equation (24) with ik, and 



using the incompressibility condition ( 23 ) , gives 

k- (Me) 



fj = 2v Q Xie\\ +i- 



Ikl 2 



Inserting this into Equation (|24|) and defining Mj 
%(k) 



II(k)M, where 



kikj 

5ij - 7^4 



is the orthogonal projector of k, we obtain 
are = — Mj_ e. 

The eigenvalue spectrum of the matrix Mi is given by 



<r(k) g ^o, r |k| 2 + r 2 |k| 



k 2 

2A^ 



i\ v k x 



(26) 



(27) 



(28) 



(29) 



The zero eigenvalues correspond to the Goldstone modes. The non-zero eigenvalues have 
eigenvectors (—k y , k x ), implying that, for T < 0, there will be a range of exponentially 
growing modes in the direction perpendicular to k. 



Equations (18) and (29) predict that, when A < and T < 0, isotropic and 



polar fixed points become simultaneously unstable, thereby signaling the existence 
of spatially inhomogeneous dynamic attractors. More generally, within the class of 



standard PDEs, the two T-terms in Equation (14) appear to provide the simplest 'linear 
way' of obtaining a v-only theory that exhibits non-trivial stationary dynamics. In 
principle, one could also try to model instabilities by combining odd or fractional 
powers of |k| in Equations (18) and (29); this would be analogous to replacing the 



quartic Landau potential by a more general function of |v|. However, when considering 
eigenvalue spectra based on odd or non-integer powers of |k|, the underlying dynamical 
equations in position space would become fractional PDEs. Such fractional models could 
potentially be useful for describing active suspensions with long-range or other types of 
more complex interactions, but their analysis goes far beyond the scope of this paper. 



3.3. Numerical results in 2D 



We simulated the vector model, defined by Equations m\ and (14), in two space- 



dimensions using an algorithm similar to that described in Section 2^ The primary 
difference compared with the simulations for the scalar model is an additional pressure 
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model parameter 


rescaled dimensionless parameter 


A 


At u 


C 


1 


r 


To t u L 2 


r 2 


1 


s 


S 



Table 2. The right column shows the rescaled dimensionless parameters used in the 
simulations of the vector model from Equation ( 14 ) . The unit time is defined by the 
damping time-scale t u = L /T%, where L is the length of the 2D simulation box, and 
the order-parameter is measured in units of v u — 1/ y/t u C. 



correction subroutine that ensures the incompressibility of the flow (see Reference (23] 
for details). Table [2] summarizes characteristic units and rescaled parameters as adopted 
in the computations. All simulations were initiated with random initial conditions, and 
the typical time discretization was At = 0.1 (in characteristic time units t u ). 

Kinematic transitions We first study how a decrease of the 'viscosity' parameter T 
affects the stationary dynamics for mono-stable and mexican-hat (polar-ordering) 
potentials. To this end, simulations were performed with fixed potential functions 
at three different values of the pusher/puller parameter S, while varying To between 
successive runs. Changes in the quasi-stationary dynamics are quantified by measuring 
the space-time averaged variance a 2 = (v 2 ) — (v) 2 , shown in Figure [3j Similar to the 
scalar model, the first transition from an isotropic state with v = to a non-trivial 
stationary dynamics with a 2 > is found to occur at T « — (27i) 2 T 2 / L 2 . 




01234 5 01234 

-r„ -r 



Figure 3. First few kinematic transitions in the quasi-stationary dynamics of the 
vector model for (a) mono-stable and (b) mexican-hat potentials at different values 
of S. Transitions are indicated by sudden changes of order-parameter variance 
<7 2 — (v 2 ) — (v) 2 . The space-time averages (•) were measured from 3000 successive 
snapshots (At = 0.1) after an initial relaxation period of 200 time units t u . 
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Figure 4. Simulation snapshots (t — 400) of streamlines and vorticity (scaled by the 
maximum value w m ) for the sombrero (polar-ordering) potential and random initial 
conditions. (a,b) Stationary vortex lattices obtained in the special case S = 1, 
corresponding to a vanishing convective derivative in Equation (14 1. (c,d) Strong 
convection (\S\ 3> 1) leads to the formation of dynamic patterns. 



To illustrate kinematic changes in the flow dynamics in more detail, we show quasi- 
stationary snapshots from simulations with a mexican-hat potential for different values 
(To, £7) in Figure |4| In the special case £7 = 1, corresponding to a vanishing convective 



derivative in Equation (14), stationary cubic vortex lattices form, with an increasing 
number of vortices as r is decreased (Figure |4^,b). By contrast, for £7^1, nonlinear 
convective effects cause distortions in the vortex lattices. As a consequence, the 
dynamical system no longer approaches a time-independent stationary state but instead 
exhibits a complex non-equilibrium dynamics (Figure j4j?,d). For pushers (£7 < 0), the 
resulting turbulent flow patterns look remarkably similar to those observed in dense 
suspensions of Bacillus subtilis [7| l8l l23| 154]. 

Flow Spectra To resolve the structure of the flow fields obtained from Equations (|8 



and (14), we calculated the energy spectrum E(k), formally defined by 

POO 

(v 2 ) = 2 E(k)dk, (30) 







where k = |k|. By virtue of the Wiener-Khinchine theorem [55], E(k) can be estimated 
by Fourier-transformation of the equal-time two-point velocity correlation function, 
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Figure 5. Flow energy spectra in arbitrary units for (a) mono-stable and (b) mexican- 
hat potentials and different values of S. Spectra are scaled to have the same maximal 
value. In the strongly non-linear regime, corresponding to |5| 3> 1, the spectra of both 
pullers (S > 0) and pushers (S < 0) seem to approach 'universal' limit functions, the 
exact shape of which depends on the type of the potential. Simulation parameters: 
grid size 256 x 256, r = —1 in all simulations. 



yielding in d dimensions 

E d (k) ~ k d ~ l J d d R e~ ikR (v(t, r) • v(t, r + R)). (31) 

Traditionally, spectral flow analysis has been an important tool in the investigation 
of classical turbulence phenomena [53] . Flow spectra for our numerical data are 
summarized in Figure [5j The critical case 5 = 1 provides a useful reference point, since 
in this case the stationary dynamics becomes static, as illustrated in Figure |4^,,b for the 
sombrero potential. Accordingly, the spectra for 5 = 1 exhibit a sharp peak that reflects 
the typical vortex size in the stationary state (see green curves in Figure [5j*,b). In the 
case of the mono-stable potential, changing the pusher/puller parameter to values 5^1 
affects the asymptotic slopes of the spectrum but leaves the position of the maximum 
practically constant (Figure [5^l). By contrast, for the polar-ordering potential, both 
slope and position of the maximum change as 5 is decreased or increased from unity 
(Figure [HJd) . The large- 15| pusher-spectra for the mexican-hat potential agree well with 
those measured in dense quasi-2D B. subtilis [23] suspensions. 

Bionematic jets To illustrate qualitatively how the velocity potentials affect the 
stationary dynamics - and, hence, the spectral flow properties - we present in 
Figure [6] snapshots from computations with an intermediate-size simulation volume 
and Tq — (27r) 2 T2/L 2 . For the mono-stable potential, we observe fairly homogeneous 
vortex structures (Figure |6^,b) in agreement with the relatively sharp spectral peaks 
in Figure [5^. By contrast, in the case of the mexican-hat potential (Figure |6]3,d), 
extended jet-like regions form that look very similar to the 'zooming bionematic' phases 
in dense B. subtilis suspensions [71 151 153] . 
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Figure 6. Snapshots of stationary flow and vorticity patterns, scaled by the maximum 
value w m , for (a,b) mono-stable and (c,d) mexican-hat potentials. The pusher flow 
field in (c) agrees qualitatively with experimentally observed flow fields for dense 
B. subtilis [3 |H1 H31 [M] suspensions. 



Nonlinear limit & 'universality 7 Interestingly, our numerical results suggest that the 
spectra approach universal functional forms in the nonlinear limit |Ao| = \S — 1| ^ 1, 
even though the exact asymptotic behavior depends on the type of the potential. 
Moreover, for pushers with A 1, the spectra obtained from Equation (14) agree 
well with those measured for dense quasi-2D B. subtilis suspensions [23], and a very 
similar small-A; scaling was also observed in 2D particle simulations of self-propelled 
rods [23]. These observations hint at some 'universality' in the spectral properties of 
dense active particle systems. Future investigations of this question will require larger 
simulations as well as systematic experimental investigations of larger systems, which 
allow to extract asymptotic scaling laws and other spectral details with higher accuracy. 
Generally, in the strongly nonlinear regime |Ao| = |1 — S\ ^> 1, the vector model predicts 
phenomenologically similar behavior for both pullers (S > 0) and pushers (S < 0), 
see Figure [6j However, while values of Ao 3> 1 appear to have been realized in dense 
suspensions of pusher-type B. subtilis bacteria [HI [221 03], it seems unclear at present 
whether puller (e.g. algal) suspensions with S ^> 1 can be achieved in experiments. 
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4. Conclusions and future challenges 

Identifying 'universal' features of active fluids is one of the key challenges en route to 
a more systematic physical classification of biological matter. Here, we have argued 
that, in many cases, phenomenological aspects of dynamical non-equilibrium phases 
can be naturally described by minimal models that focus on a limited number of 
experimentally accessible order-parameters and are based on higher-than-second-order 
spatial derivatives. If true, this would imply that a variety of active fluids with 
different microscopic constituents can exhibit very similar long-wavelength behavior. 
More generally, adopting the view that kinematic transitions in living fluids reflect 
qualitative changes in small-wavenumber expansions, and thus may be interpreted as 
Landau-type transitions in Fourier space, can help to catalog non-equilibrium systems 
according to their asymptotic behavior at long wavelengths, similar to the classification 
of phase transitions in equilibrium thermodynamics. 

Simplified continuum descriptions of active fluids remain practically relevant 
because many of the recently proposed multi-order-parameter theories feature a 
large number of unknown transport coefficients and, therefore, cannot be tested 
in detail with present and next generation data. While our above analysis has 
focussed on the most basic scalar and vector models, the approach can be easily 
extended to higher tensorial order-parameter fields (e.g., Q-tensor descriptions [56] 
of active nematics). To test whether microscopically different active fluids do indeed 
share universal hydrodynamic long-wavelength characteristics will require combined 
analytical, numerical and experimental efforts. First attempts to apply the above ideas 
to very dense bacterial suspensions give promising results [23] but further quantitative 
studies will be needed to decide whether fourth and higher-order PDEs are capable 
of providing a sufficiently accurate phenomenological description of experimentally 
observed active phases. To expedite quantitative comparisons with experiments, it 
would be desirable to develop alternative computation schemes that will allow for 
faster simulations of higher-order PDE models. Promising candidates could be suitably 
adapted Lattice-Boltzmann algorithms [57] that implement negative 'viscosities' [58] . 

Finally, the theoretical analysis of higher-order PDEs poses a number of 
future mathematical challenges, one of which being their derivation from underlying 
microscopic, multi-component or kinetic models through systematic projection methods. 
Another particularly interesting question relates to the formulation of consistent 
boundary conditions at interfaces. This problem, which was circumvented in our 
simulations by considering toroidal domains, is also encountered in other higher-order 
structure formation and transport theories [291 EHJ EE] • To identify physically reasonable 
boundary conditions for active fluids at solid interfaces is not only relevant from 
a mathematical perspective but also from the experimental point of view, as they 
will directly affect predictions for effective shear viscosities |59j and other measurable 
quantities. 
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